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Abstract 

We present an extension of a simple automaton model to incorporate non-local 
interactions extending over a spatial range in lattice gases. From the viewpoint of 
Statistical Mechanics, the lattice gas with interaction range may serve as a prototype 
for non-ideal gas behavior. From the density fluctuations correlation function, we 
obtain a quantity which is identified as a potential of mean force. Equilibrium and 
transport properties are computed theoretically and by numerical simulations to 
establish the validity of the model at macroscopic scale. 

KEY WORDS: Lattice gas automata; interaction potential; fluctuations correlation 

function; spinodal decomposition. 
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1 Lattice gas with non-local interactions 

Standard lattice gas automata (LGA) evolve according to an iterated sequence of mass- 
and momentum-preserved local collisions followed by propagation. Non-local interac- 
tions can be incorporated in the LGA dynamics via long-distance momentum transfer 
simulating attraction and/or repulsion between particles H, ||, ||, In local collisions, 
momentum redistribution is a node-located process with local conservation of mass and 
momentum. In non-local interactions (NLI), momentum is exchanged between two par- 
ticles residing on nodes separated by a (fixed or variable) distance r: mass is conserved 
locally, momentum is conserved globally. At the macroscopic level, the main feature 
exhibited by LGA models with NLI's is a "liquid- gas" -type phase separation with bub- 
ble and drop formation [I]. From the statistical mechanical viewpoint, LGA's with NLI's 
form an interesting class of models in that — in contrast to standard collision-propagation 
models — they include an elementary process which is essential for "non-ideal" behavior. 

The dynamics of LGA virtual particles is not governed by Newton's equation of motion 
and the concepts of force and potential cannot be used in the sense of classical mechanics. 
Moreover in real fluids each particle is subjected a priori to the force field of all particles 
(whose effect is quantified by the potential of mean force) whereas in discrete lattice gases 
each particle interacts non-locally with at most one other particle at a time. So stricto 
sensu the usual concept of intermolecular potential does not apply to lattice gases. 

In the LGA model with NLI's proposed by Tribel and Boon [Q, the idea of an interac- 
tion range was introduced by governing the interaction distance according to a probability 
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distribution — namely a power law (oc r~ M ) — wherefrom a distance r is drawn for each 
particle at every time step. Here we show that for sufficiently long times and large number 
of particles, the implementation of a probability distribution of interaction distances has 
a resulting effect similar to the effect of an interaction potential. We first describe the 
model in section |2|. Then in section § we compute the density fluctuation correlations §| 
wherefrom a quantity is obtained which can be identified as a potential of mean force. 
Sections f| and [| present an analysis of the equilibrium and transport properties. We 
conclude with some comments. 



2 Interaction range model 

The automaton resides on a two-dimensional triangular lattice and uses for propagation 
and local collisions the rules of the FHP-III model |4j with periodic boundary conditions. 
Non-local interactions can take place between two particles when nodes separated by some 
distance r exhibit favorable configurations as illustrated in Fig.l. The interaction modifies 
the orientation of the velocity vectors from a diverging configuration to a converging 
configuration to simulate attractive forces and vice-versa for repulsive forces. At each 
time step, the algorithmic procedure must realize a pairing of particles separated by a 
distance r drawn from a probability distribution p(r). It is clear that a parallel algorithm 
can hardly be efficient here. Therefore we use a sequential algorithm which proceeds as 
follows: 

(i). at each time step, a direction is arbitrarily chosen along any of the lattice axes and 
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all interactions will be along that direction during that time step; 

(ii) . a particle, say at node A, is (sequentially) selected and accepted if its state has not 

been modified by a previous interaction in the sequential procedure; 

(iii) . a distance r is drawn from the distribution p(r) and a pointer is set at nodes F and 

B, located respectively at a distance r forward and backward from A; 

(iv) . if one of the configurations 11 B A" or "AF" is compatible for interaction (see Fig.l), 

the configuration is modified accordingly and the procedure keeps track of the mod- 
ification for the duration of the sequence (each particle can undergo no more than 
one interaction per time step). 

As a result the effective probability that an interaction occurs in the simulation differs 
from the theoretical p(r). The details of the computation are given in the Appendix; here 
we merely quote the final result which expresses the effective probability q(r) in terms of 
Pf(j') and psij'), denoting respectively the forward and backward probabilities with the 
imposed analytical form (e.g. p(r) oc r~ M ): 

q(r) = p F (r) + p B (r) - p F (r)p B (r) (1) 



with 



and 



p F (r)=p(r) J] [1-«2Pf(*)] (2) 

£=r+l 



{fmax ^ ^ r 1 

n [i-K 2PF (£)}\ ni 1 ( 3 ) 



where r max is the cutoff distance in the distribution p(r) (see [§]) and k 2 = /(l — /), 
with / the particle density per channel 1 . Besides the fact that here interaction distances 
are distributed over an interaction range, an important difference with the fixed-distance 
model (y|, [§]) is that in the present case, each particle belongs to many pairs of possible 
interactions. Note also that for similar reasons, there is a bias in the effective distribution 
toward large interaction distances. Indeed when drawing a low value of r from p(r) in 
the sequential procedure, there is a greater chance that the second particle of the pair 
be already involved in a previous pairing. As a result the effective probability for long 
distance interactions is larger than predicted by the pre-set distribution p(r). 

The question now arises as to define a quantity which can be identified as an interaction 
potential in a discrete lattice gas. We propose the following heuristic argument. We 
evaluate the rate of momentum exchange caused by the non-local interaction 

F(r) = l4q(r) (4) 

where 7 is a numerical factor whose value corresponds to the average amount of momen- 
tum transfer (7 = 4/3 and 7 = 1 for the models shown in Figs, l.a and l.b respectively). 
Interpreting F(r) in (f|) as a force, we define the "pair potential" as the discrete analogue 
of the potential in continuum mechanics: 

u(r) = -Kf±F(£) 
e=i 

= -E7?W 

t 

1 Note that if desired operationally, Eqs. (0)-(§) can be inverted numerically to obtain a function p(r) 
such that the effective distribution q(r) be of given analytical form. 
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= -iHr) (5) 

where T{r) is the repartition function corresponding to the distribution q(r). Then using 
Eqs. (P])-(H), u(r) is well-defined once p(r) is fixed. For instance, if we use the power-law 
distribution p(r) oc r~ M such that the interactions are repulsive for r = 1 and attractive for 
r = 2, . . . , r max , u(r) exhibits a form compatible with the expected typical pair interaction 
potential, as shown in Fig. 2. 



3 Density fluctuation correlations 

The next question is the influence of the non-local interactions on the density fluctuation 
correlations in the lattice gas, which is most conveniently measured by the static structure 
factor || defined by 

^)4lE fe !(Mm(M), (6) 
t=l i,j 

where p is the density per node, and 

6m(k,t) =$>- lk - x Kx,c i; t) - /] (7) 

X 

is the fluctuation of the channel occupation number rij (i = . . . b). In the ideal lattice gas 
(whose dynamics is governed by propagation-collision rules) there are no static density 
correlations and the static structure factor is a constant ||: 

S°(*) = (l-/)(1 -<$(*)). (8) 

By analogy with the statistical mechanical theory of continuous fluids , we write 

S(k) 



S°(k) 
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1 + fh(k) (9) 



where h(k) is the Fourier transform of the pair correlation function [g(r) — 1] and is 
therefore related to the potential of mean force <p(r) since g(r) = exp[— /3<p(r)] (here 
(3 is an arbitrary constant). So by measuring the density fluctuation correlations in 
lattice gas simulations, we can extract a function 0(r) from the measured static structure 
factor. The results are shown in Fig.3: both the radial distribution function g(r) and the 
potential function <p(r) are reminiscent of those obtained in real fluid measurements @. 
The connection between the potential of mean force <j>(r) and the interaction potential 
u{r) discussed in section || remains to be clarified. 

Another interesting feature is worth mentioning. Consider the LGA is in the appropri- 
ate density range for spinodal decomposition (see section [|). Then one could effectively 
"quench" the system by increasing the interaction range. By measuring S(k) at suc- 
cessively increasing values of r max we find that S(k) increases dramatically at low k. 
Following the lines of heuristic reasoning and anticipating a result of section we infer 
that 

g(fc-»0)~ 1-/ (10) 
l-jK 3 {r) q 

where the denominator follows from the expression for the compressibility (see section 
f|). Here (r) q is the expectation of r computed with the distribution q(r). Since (r) q 
increases with r max , S(k — * 0) grows accordingly as expected when the phase transition 
is approached. Further analysis will be presented in a forthcoming paper. 
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4 Equilibrium properties 

The pressure at global equilibrium is given by (see ||) 

p4EEK( x ) + m .(4 (ii) 

where V is the number of nodes of the automaton universe, J2xgc Si(^i( x )) is the mo- 
mentum transport due to propagation and Sxe£ Si( m i( x )) is the momentum flux due to 
NLI's. Then the hydrostatic pressure can be evaluated as follows: 

• the convective momentum flux is the total momentum carried by moving particles 
in the fluid. On each node in the FHP-III model, there are 6 channels with velocity 
1 and one zero-velocity channel, so that 

££M*)> = 4a (12) 

xe£ i 

where p is the average density per node; 

• the non-local momentum flux is caused by NLI's. The value of this flux is clearly 
given by 

££fa(x)) = 3V7K2£Hz(r), (13) 

xe£ i r 

since on the average each NLI causes a momentum flux of value jr and on a given 
node 3 particles may independently be involved in an interaction. 

Consequently 

P=V-\l4(r) q (14) 
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with / = f and (r) q = J2r rq(r). The first term on the r.h.s. of (|H| ) is the kinetic pressure 
of the ideal lattice gas and the second term as given in this mean-field evaluation depends 
only on the first moment of the distance distribution q{r). Note that for fixed-distance 



interaction models, q{r) = 5(r — £) and ( |l4l) becomes for the model of Fig.l.a (with 7=0) 



p = 3f- 2£k 2 2 (15) 



as given in and [§. 

From ( |14"|) it follows that the compressibility is given by 

I dp 
p op 

1 

fp(l-7(r),K 3 ) 


= 1 ~, v , (16) 
1 -7(r) ? K 3 

where k 3 = /(l — /)(1 — 2/) and x° is the compressibility of the ideal gas. By us- 
ing S{k — > ) = p/3~ 1 Xth |1 and the thermodynamic pressure of the ideal gas p t h = 
— 6/3 _1 ln(l — /) which yields the compressibility Xa = P~ l Pi^ ~ f) (with = Cq = 3/7 
for the FHP-III model) we obtain Eq.fllJ). 

The compressibility equation ( |i4|) yields the square of the sound velocity c s 

^ = ^ = cg[l-7«8<r>J. (17) 

Eq. ( |TTD is valid as long as |^ > 0; when |^ < 0, the density fluctuations show an 
explosive behavior and the system separates into two phases (i.e. for (r) q > 7.79). In Fig. 4 
we show the results of measurements of the pressure as given by Eq. (0), compared to 
the theoretical prediction fli"4|) ; in Fig.5, the theoretical sound velocity (|17D is compared to 



the simulation results (for the experimental method, see e.g. 0). Experimental evidence 
of spinodal decomposition is given in Fig. 6 where the evolution of the density distribution 
shows how phase separation takes place in the automaton. 

5 Transport coefficients 

5.1 Microdynamical and lattice Boltzmann equations 

The evolution of the automaton is obtained by applying successively the non-local inter- 
action routine, the collision routine, and the propagation routine. This computational 
procedure is the operational realization of the microscopic dynamics of the automaton 
whose mathematical formulation is given by the microdynamical equation 

n(x + c», Ci] t + l)=C l {X [n(x; *)]} (18) 

where n(x, Cjji) is the Boolean occupation variable of channel i at node x at time t. C 
and X are the local collision and non-local interaction operators respectively. The explicit 
expression of the non-local operator X reads 

6 r i=-l 

with (for the model of Fig.l.a) 2 

^i n ( x ; t) = [n;(x; t)n i+3 (x; £)] [rij(x + rcf, £)n i+3 (x + rcf, t)] 
2 Channel indices are numbered counter-clockwise from to 5 for moving particles and 6 for the rest 

particle, and indices i and i + j are taken modulo 6. 
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- [n 4 (x; t)n i+3 (x; £)] [rTi(x - rc^; t)n i+3 (x - rc^; t)] 
2^ ±1 n(x; t) = [rTi(x; t)n i=Fl (x; i)] [n^x + rcy-i; t)n i=F i(x + rc i± i; *)] 

- [rii(x; t)7T iT i(x; £)] [rii(x - rc i± i)n i=F i(x - rc i±1 ; £)] , (20) 

where 7T = 1 — n. 

Taking the average of Eq. ([T^) over a non-equilibrium ensemble, and making the 
molecular chaos assumption, one obtains the lattice Boltzmann equation 

/(x + a, a; t + 1) = d {1 [/(x; *)]} , (21) 

where /(x, Cjji) = (n(x, cj; £)) is the singlet distribution function of channel i at node x 
at time i. In this equation, the operators C and X act on the distribution function / (not 
on the Boolean variables rij). 

5.2 Linearized lattice Boltzmann equation 

Considering small deviations from local equilibrium 

f{x,c i] t) = f + 5f(x,* ] t) (22) 

the lattice Boltzmann equation ( PT{ ) may be linearized for the perturbation Sf. Denoting 
by Q the usual linearized collision operator and by A the linearized NLI operator 

5f(x,c i ;t) = (l + A) ij 5f(x,c;t), (23) 

the linearized lattice Boltzmann equation reads 

5f(x + a; Cij t + l) = (l + n) -(l + A) jk 5f(x, c k ; t). (24) 
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We develop the perturbations 5f in Fourier modes 



l" k 



(25) 



and rewrite Eq.([24j) in Fourier space to obtain the eigenvalue equation for the automaton 



(26) 



This equation is formally identical to the eigenvalue equation for the fixed-distance model 
0, but the operator A is now a linear combination of the corresponding fixed-distance 
operators. The three slow (hydrodynamic) modes of interest are the shear mode, noted 
ip u , whose eigenvalue corresponds to the kinematic viscosity v, and the two sound modes 
ip a= ± related to the sound velocity c s and damping coefficient T. 

5.3 Transport coefficients 

As mentioned, the operator A is a linear combination of fixed-distance NLI operators 



K = E <?( r ) A 



*7* 



(27) 



where A*J denotes the NLI operators at distance r. So the results for the fixed-distance 
interaction operators can be extended to the present model (see for details). The main 
difference is that the distance r is now replaced by its mean value (r) 9 , and r 2 by the 
variance (r 2 ) q = J2 r r 2 q(r). Expanding the eigenvalues in powers of k 



s M (k) = (ik)z<f> + {ikfz^ + 



(28) 
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we obtain 



l x) = ±c s , 4 1} = o, 

(29) 

(2) 



4 2) = r = ±(i/ + c), z r = v. 



with 



v = "0 U - ~(r) 9 K 3 ) (1 - (r) g K 3 ) + J^( r )i K 3 U - ~(r) q K^\ + ^(r 2 ) q K 2 , 

C = Co{l-l(r) q ^-^(r) q ^ + \(r 2 ) qK2 . (30) 

vq et Co are the kinematic and bulk viscosities of the standard FHP-III model 

1/1 1 



U ° 14 \u„ 2 

* ■ MH- 

with 

^i/ = ^2(7-8k 2 ), 

cu c = 7« 2 (l-2« 2 ). (32) 

Excellent agreement is obtained between the simulation data and the theoretical results 
as shown in Fig. 7. 



6 Comments 

The question was raised by Gerits et al. || that models with non-local fixed interaction 
distance lack detailed balance and that therefore their equilibrium distribution is not 
known. In the model presented here the constaint is weaker in that interaction distances 
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are distributed and the non-local interactions are probabilistic. Yet the condition for semi- 
detailed balance was not taken into account. The mean field theory is found to predict 
correctly macroscopic equilibrium and transport properties. Fluctuation correlations were 
measured and the static structure factor was used to extract the LGA analog of a potential 
of mean force. The statistical mechanical theory for the static and dynamic structure 
factors will be discussed in a forthcoming paper. 

A Evaluation of q(r) 

Here we show how to compute the effective automaton distance distribution q(r) from 
a given probability distribution p(r). The algorithmic procedure considers each node 
sequentially. Define the node currently examined as the "center node" A and define the 
"forward node" F and "backward node" B located along the direction of interaction at 
a distance r on each side of A. Now each particle on A may interact with at most one 
particle located either on F or on B. However since the algorithm is sequential, the 
forward and backward probabilities are different: a backward interaction is possible only 
if the particle on A has not interacted before, while a forward interaction is independent 
of previous interactions involving node A. 

• Forward interaction probability : Pf(j')- Suppose that configurations on A and F 
are favorable (see Fig. 1). Then the only additional condition is that the particle 
on F cannot have been involved previously in an interaction from a distance larger 
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than r. This "non-event" has the probability 



(33) 



t=T + l 



Consequently the equation for pp is given by 



p F (r) = p(r) Jl [1 - k 2 Pf{£)\ ■ 



(34) 



l=r+l 



• Backward interaction probability : ps(r). We must consider that (i) the particle 
on A has not been involved in a previous (forward) interaction, (ii) the particle 
on B has not been paired successfully with another particle when B was a center 
node (forward interaction), and (iii) the particle on B has not been paired success- 
fully with a center node located between B and A (backward interaction). The 
corresponding probabilities are: 

(i). no interaction with A forward node: 




(35) 




(iii). no interaction with B as backward node: 



r-l 



n 




(36) 



The backward probability is therefore 




(37) 
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As a result q(r) is expressed as a combination of p B and p F , given that forward and 
backward interactions are mutually exclusive: 

q(r)=p F (r)+p B (r)[l-p F (r)]. (38) 

In the case of fixed-distance interactions (at a distance £ ), the above considerations do 
not apply and the distribution reduces to 

q(r) = 5 rA) . (39) 
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Figure Captions. 

Figure 1: Interaction configurations: illustration of configuration changes through non- 
local interactions between pairs of particles on nodes at distance r from each other. Dotted 
(full) arrows indicate channel occupation before (after) interaction: [Qj,X k ] — > [Q k ,Xj], 
X = R, S,T where the channel indices (j, k) are given modulo 6 for % — 0, . . . , 5. Momen- 
tum exchange through interactions is two or four units for (a) configurations whereas all 
(b) interactions exchange two momentum units. 

Figure 2: Pair interaction potential: u(r) for a power-law distribution p(r) oc r~ M for 
1 < r < r max , with fj, = 1.2 and r max = 10. Potential units are arbitrary. 

Figure 3: (a) Radial distribution function g(r). (b) Potential function <3>(r) = \n(g(r)). 
Probability distribution p(r) oc for 1 < r < r max ; r max — 6, /i — (circles), r max = 8, 
/i — (squares), r max — 10, /i — 1 (diamonds). Lines are guides to the eye. Mean channel 
density / = 0.1; lattice size 512x512; g{r) measured over 500 time steps. 

Figure 4: Pressure versus density, (a): fixed-distance interactions; (b): distributed- 
distance interactions, with "flat" distribution p(r) oc r^ ax . Symbols are experimental 
data, curves are theoretical predictions. Fixed distances r and cutoff distances r max are 
as indicated. Lattice size 512x512; each point obtained by averaging over 300 time steps. 

Figure 5: Sound velocity versus density. Circles, squares, diamonds (resp. full, dotted, 
and dashed lines) correspond to fixed- distance interactions; triangles (long-dashed line) 
correspond to "flat" distributed-distance interactions. Symbols are experimental data, 
curves are theoretical predictions. Fixed distances r and cutoff distances r max are as 
indicated. Lattice size 512x512; each point is obtained by averaging over 10 runs of 2500 
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time steps. 

Figure 6: Evolution of the density distribution, measured every 100 time steps for a 
total duration of 700 time steps. The evolution shows horizontal separation of density 
peaks characteristic of spinodal decomposition, as opposed to vertical growth of peaks in 
nucleation and growth processes. Lattice size 512x512; /i — 0; r max = 20; density values 
are coarse-grained by averaging over each node and its six nearest neighbors. 

Figure 7: (a) Kinematic viscosity v versus density; (b) Sound damping coefficient T 
versus density. Symbols are experimental data, curves are theoretical predictions, condi- 
tions are as indicated. Lattice size 512x512; each point obtained by averaging over 100 
runs of 2500 time steps. 
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